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Q 2 evolution equations are important not only for describing hadron reactions in accel- 
erator experiments but also for investigating ultrahigh-energy cosmic rays. The standard 
ones are called DGLAP evolution equations, which are integrodifferential equations. There 
are methods for solving the Q 2 evolution equations for parton-distribution and fragmenta- 
tion functions. Because the equations cannot be solved analytically, various methods have 
been developed for the numerical solution. We compare brute-force, Laguerre-polynomial, 
and Mellin-transformation methods particularly by focusing on the numerical accuracy 
and computational efficiency An efficient solution could be used, for example, in the 
studies of a top-down scenario for the ultrahigh-energy cosmic rays. 

o 

O ; 1. Introduction 

Q_i! High-energy hadron reactions are described in terms of parton-distribution functions 

(PDFs) and fragmentation functions (FFs). There are parametrizations for the PDFs [ 
d and FFs [12] . Using these funct ions, cross sections of high-energy hadron reactions are 
evaluated. Precise calculations of these cross sections are important for finding any new 
physics beyond the current theoretical framework. 

The PDFs depend on two kinematical variables x and Q 2 . They are defined by Q 2 = —q 2 
and x = Q 2 /(2p ■ q) in lepton scattering with the momentum transfer q and the hadron 
momentum p. The FFs depend on Q 2 and another variable x = 2Eh/\fs, where Eh is the 
hadron energy and ^/s is the center-of-mass energy. Their Q 2 dependence is called scaling 
violation, which is calculated by the DGLAP (Dokshitzer-Gribov-Lipatov-Altarelli-Parisi) 
evolution equations [Ej in the perturbative QCD region. 

The Q 2 evolution equations are frequently used in describing high-energy hadron reac- 
tions. Because the PDFs and FFs vary significantly in the current accelerator-reaction 
range, Q 2 =l GeV 2 to 10 5 GeV 2 , the Q 2 dependence should be calculated accurately. Fur- 
thermore, it is known that high-energy cosmic rays have energies much more than the 
TeV scale. Analytical forms of current PDFs and FFs are supplied typically in the GeV 
region, so that they have to be evolved to the scale which could be more than TeV in 
order to use them for investigating the cosmic rays [IHE]. 

A useful evolution code was developed in Ref. [ H] for the cosmic ray studies. The 
Laguerre-polynomial method was used for solving the evolution equations. Splitting func- 
tions, PDFs, and FFs are expanded in terms of the Laguerre polynomials, and then the 
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evolution is described by a simple summation of their expansion coefficients. In fact, the 
method is very efficient for solving the equations in comparison, for example, with a direct 
integration method [EH El- However, because the Laguerre polynomials L n {— lnx) go to 
infinity in the limit x — > 0, it could have an accuracy problem in the small x region where 
high-energy reactions are sensitive. In this paper, we discuss evolution results by the 
Laguerre method [HJ|S] in comparison with the ones by other solution methods, "brute- 
force" [IH] and Mellin-transformation [ 9 methods. In particular, evolution accuracy and 
computation time are compared. It is the purpose of this paper to clarify the advan- 
tages and disadvantages of these numerical solution methods for a better description of 
high-energy hadron reactions including the high-energy cosmic rays. In particular, the 
FF evolution could be used for studying a top-down scenario in order to determine the 
origin of ultrahigh-energy cosmic rays, namely from a decay of a superheavy particle [|5|. 

This paper consists of the following. The DGLAP evolution equations are introduced 
in Sec. El and numerical solution methods are explained in Sec. El Evolution results and 
their comparisons are discussed in Sec.EJ The results are summarized in Sec. 03 

2. Q 2 evolution equations 

From the cross-section measurements of high-energy lepton-hadron, hadron- hadron, 
and lepton-annihilation reactions, the PDFs and FFs are extracted. The PDFs and FFs 
are expressed in terms of the two kinematical variables x and Q 2 . A PDF or a FF is 
expressed f(x, Q 2 ) in the following. We investigate the standard Q 2 evolution equations, 
which are called the DGLAP evolution equations [E]. The flavor nonsinglet equation is 
written as 



d\nQ 2 



L s (^Q 2 ) = ^r 1 [ d -^P NS {x/y)f NS {y,Q 2 ), (i) 

2vr J x y 



where f NS (x,Q ) is a nonsinglet (NS) function, P NS (x) is a nonsinglet splitting function, 
and a s (Q 2 ) is the running coupling constant. The splitting functions for parton distribu- 
tions and fragmentation functions are identical in the leading order (LO) of a s ; however, 
they differ if higher order corrections are included [CSj. In order to make the evolution 
equation slightly simpler, the variable t is used instead of Q 2 : 



2 , 



a s (Q 2 ) 
a s {Ql) 



(2) 



where the running coupling constant in the leading order (LO) is given by a s (Q 2 ) = 
47r/[/3o ln(Q 2 /A 2 )] with the QCD scale parameter A. The constant (3q is expressed as 
fa = HC G /3 - AT R N f /3 with C G = N c and T R = 1/2. Here, N c is the number of color 
(N c =3) and Nf is the number of flavor. The Qq in Eq. (J2J) indicates the initial Q 2 where 
the evolved function is provided. Using this variable t in Eq. (fT|). we obtain 

^firs( x > t ) = J ^ p Ns{ x /y) fN S (y,t), ( 3 ) 

for the nonsinglet evolution. 
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In the flavor singlet case, the evolution is described by two coupled integrodifferential 
equations: 



d_ 
dt 



f(x,t) 



1 dy 

y 



P(x/y) f(y,t), 



where the matrices f and P are defined by 



f{x,t) 



fs( x ^) 
fgfat) 



x 



Pqq \X) 



2NfP ij {x) 

Pgg ( x ) 



(4) 



(5) 



The indices i and j indicate ij = qg for the PDFs and ij = gq for the FFs. The functions 
Pqq, Pqg, Pgq, and P gg are splitting functions. The function P^ determines the probability 
of the splitting process that a parton j with the momentum fraction y splits into a parton 
i with the momentum fraction x and another parton and the j-parton momentum is 
reduced by the fraction z. In the LO, the splitting functions are expressed as 



Pqq( x ) —Pns( X ) 



Ct, 



1 + x 2 3 r/1 , 



Pqq(x) =T R [X 2 + (1-X) 2 ], 

l + (l-x) 2 



Pgq(x) =C F 



P gg (x) =2C G 



X 



(6) 



x 



1-x) 



+ 



1 — X 



X 



+ x(l-x)+ [ — --^f^ ) 6(1 -x) 



12 3 a 



G 



where Cp is given by Cf = (N 2 — l)/(2iV c ) and 1/(1 — x) + is defined by dxg(x)/(l — 
x) + = Jq 1 dx [g{x] — <?(1)]/(1 — x) with an arbitrary function g(x). 

We need to solve the nonsinglet and singlet evolution equations in Eqs. © and (jlj. 
These are not simple integrodifferential equations, so that an efficient numerical method 
should be investigated. In the next section, three popular numerical methods are ex- 
plained. 



3. Numerical Methods for solving Q 2 evolution equations 

There are various ways for solving the DGLAP equations. In this section, we explain 
three popular methods, brute-force, Laguerre-polynomial, and Mellin-transformation meth- 
ods. In the following subsections, only the nonsinglet evolution is explained because the 
singlet evolution can be solved in the same way. 

3.1. Brute-force method 

The simplest way is possibly to use the brute-force method [ U\ ■ It may seem to 
be too simple, but it is especially suitable for solving more complicated equations with 
higher-twist terms [E]. These equations could not be easily handled by the orthogonal- 
polynomial methods such as the Laguerre-polynomial one in Sec. 13.21 and by the Mellin- 
transformation method in Sec. 13.31 Furthermore, a computer code is so simple that the 
possibility of a program mistake is small, which means the code could be used for checking 
other numerical methods. These are the reasons why it was investigated in Ref. [|7j. 
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In the brute-force method, the two variables t and x are divided into small steps, and 
then the differentiation and integration are defined by 

df(x,t) f(x i ,t j+1 )-f(x i ,t j ) " X 



dt Atj 



/ dxf(x,t) ^2 Axkf{xk,tj), (7) 
"* k=i 



where Atj and Ax k are the steps at the positions j and k, and they are given by Atj = 
tj + i — tj and Axk = Xk — x k -i- The numbers of t and x steps are denoted N t and N x , 
respectively. Applying these equations to Eq. (J3J), we write the nonsinglet evolution from 
tj to tj+i as 

Nx Ax k 

q NS (xi,t j+1 ) = q NS (xi,tj) + Atj —— p gq(x t /x k )q NS (x k ,tj). (8) 

k=% 

If the distribution q NS is supplied at t 1 = 0, the next one q NS (x,t 2 ) can be calculated by 
the above equation. Repeating this step N t — 1 times, we obtain the final distribution at 
tx t . However, it is obvious that the step numbers N t and N x should be large enough to 
obtain an accurate evolution result. 

3.2. Laguerre polynomial method 

The evolution equations could be solved by expanding the distribution and splitting 
functions in terms of orthogonal polynomials. A popular method of this type is to use 
the Laguerre polynomials [EJIH1- They are defined in the region from to oo, so that the 
variable x should be transformed to x' by the relation x' — — In x. 

The nonsinglet evolution is discussed in the following. The evolution function E NS (x, t), 
which describes the distribution evolution from t — to t, is defined by 

f NS {x,t)= j ^E NS (x/y,t)f NS (y,t = 0). (9) 

j x y 

Then, it satisfies 

^-E NS (x,t) = ^ ^P NS {x/y)E NS {y,t). (10) 

Because this is the same integrodifferential equation as the original DGLAP equation, 
one may wonder why such a function should be introduced. There is an advantage that 
the evolution function should be the delta function at t — 0: E NS (x,t = 0) = 5(1 — x) 
because of its definition in Eq. (jHJ). It makes the following analysis simpler. The functions 
are expanded in terms of the polynomials: P NS (e~ x ) = Ylm^Ns^ 71 ^'^ anc ^ E NS (e~ x ,t) = 
Y2 n E™ (t)L n (x f ), where P™ s and E™ s (t) are the expansion coefficients. The coefficient 

F n for a function F(x) is given by F n = dxL n (x') F(x), and it could be calculated 
analytically for a simple function. If the two functions on the right-hand side of Eq. (fTTHl 
are expanded, it becomes an integration of two Laguerre polynomials. Using the formula 
Jq dy'L n (x' - y')L m (y') = L n+m (x') - L n+m+l (x') for this integration, we obtain 

j t Ks(t) = ixr-cr" 1 )^)- (ii) 

m=0 
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Because the evolution function is a delta function at t = 0, all the expansion coefficients 
are one. Therefore, this equation is easily solved to give a summation form: 



k m-1 



k=0 ' i=k 

This recursion relation is calculated with the relations B® = 1, B\ =X/i=i(-^ s — ^ivs 1 )' 
and Bq = B\ = ■ ■ ■ = B^_ x = 0. After all, the evolution is calculated by the simple 
summation: 

Nhag n 

f NS {x,t) =J2J2 \ E n~m{t) - ^n-m-l(t)] L n (-hl x)f™(t = 0). (13) 
n=0 m=0 

In this way, the integrodifferential equation becomes a simple summation of Laguerre- 
expansion coefficients, so that this method is considered to be a very efficient numerical 
method for the solution. 

3.3. Mellin transformation method 

The Mellin transformation method is one of the popular evolution methods [IHj- It 
is used because the Mellin transformation of the right-hand side of Eq. (jSJ) becomes a 
simple multiplication of two moments, namely the moments of the splitting function and 
the distribution function. The moments of the splitting functions (anomalous dimensions) 
are well known, and a simple functional form is usually assumed for the distribution at 
certain small Q 2 so as to calculate its moments easily. Then, it is straightforward to 
obtain the analytical solution in the moment space. Furthermore, the computation time 
is fairly short. These are the reasons why this method has been used as a popular method. 
For example, it is used for the x 2 analysis of experimental data for obtaining polarized 
PDFs [l!2|. whereas the brute-force method is employed in Ref. [ITS]. 

The Mellin transformation and inversion are defined by 

rl i pc+ioo 

f(s,t)= dxx s ~ l f(x,t), f( x ,t) = — dsx- s f(s,t). (14) 

Jo Jc-ioo 

Here, the upper limit of the x integration is taken one because the distribution f(x) 
vanishes in the region x > 1. The Mellin inversion is a complex integral with an arbitrary 
real constant c, which should be taken so that the integral Jq dxf(x)x c ~ 1 is absolutely 
convergent. If this transformation is used, the integrodifferential equations become very 
simple. For example, the nonsinglet evolution equation becomes 

^-f Ns (s,t)=P NS (s)f NS (s,t). (15) 
Its solution is simply given by 

f NS (s,t)=e p »s(°)t f Ns ( s , t = 0). (16) 
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Because the moments P NS (s) are well known 
quantities and the moments of the initial func- 
tion f NS (s, t — 0) could be evaluated, it is 
straightforward to calculate the evolution in 
Eq. (JTfij) in the moment space. However, 
the numerical integration is needed for the 
Mellin inversion in Eq. (J 14)) for transforming 
the moments into a corresponding x distribu- 
tion. Practically, the Mellin inversion is cal- 
culated along the integration contour in Fig. 
n Changing the complex integration variable 
s for the real one z by s = c + ze 1 ^ in Eq. (|14|) , 
we have 



(17) 

The constant c and the angle are shown in Fig. ^ Using the Gauss-Legendre quadrature 
for this integration, we obtain the evolved distribution in the x space. 

4. Comparison of evolution results 

In comparing evolution results of three methods, we take the evolved distribution of 
the brute-force (BF) method with N t =200 and A^=4000 as a standard for assessing other 
evolution accuracy. It is shown in Ref. [|7j that the evolution accuracy is better than 2% 
if N t =200 and A^=1000 are taken. This is the reason why it is taken as the standard. 
Because the details are discussed in Ref. [Ej for the evolution accuracy of the BF method, 
we discuss only the comparison with the results of the Laguerre-polynomial and Mellin- 
transformation methods. 

4.1. Parton distribution functions 

In order to show the Q 2 evolution of the PDFs, we use the MRST02 distributions [ IT4] 
which are provided analytically at Q 2 =l GeV 2 . The distributions are evolved to Q 2 =100 
GeV 2 with the MRST02 scale parameter by three evolution methods. Then, the ratio of 
the evolved distribution to the one by the brute-force method with iVt=200 and N x =A000 
is shown for finding the numerical accuracy. 

First, the evolution results of the nonsinglet distribution x(u v + d v ) are shown in 
Fig. |21 for the Laguerre method. The number of the Laguerre polynomials Ni ag is 
taken as N Lag =5, 10, 20, and 30, and each distribution ratio x(u v + d v )L aguerre / x{u v + 
d v )BF{Nt=2oo,N x =mQO) is shown. It is obvious that accurate evolution cannot be obtained 
if the number Ni ag is small in the small- and large-x regions. In particular, the ratio 
shows oscillatory behavior at small x, which results from the functional behavior of the 
Laguerre polynomials. The Laguerre polynomials L n (— In x) are shown as a function of 
x in Fig. El We find that the oscillatory functional form at small x gives rise to the 
oscillatory behavior in Fig. |21 Therefore, one should be careful in the Laguerre method 
that a large number of polynomials should be taken to obtain an accurate evolution at 
small x. Furthermore, one should be also careful in the very-large- a; region. 



Im (s) 





""\ .0^ 




/c 




/ s = c + z 



Re (s) 
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Figure 2. Evolved nonsinglet distribution ra- 
tios x(u v + d v ) Laguerre /x(u v + d v ) BF are shown 
for NLag—5, 10, 20, and 30 in the Laguerre- 
polynomial method. 



Figure 3. Laguerre polynomials L n (—hxx). 



Next, the evolution results are shown in Fig. 0]for the Mellin-transformation method. 
The Mellin inversion of Eq. (|17j) is numerically calculated by the Gauss-Legendre quadra- 
ture with the number of points Nql, which is taken as Ngl=G, 10, 20, and 50 in Fig. 
0J The integration contour of Fig. ^ is used with the constants c=l.l and 0=135° as 
suggested in Ref. [E]. In order to use the Gauss-Legendre quadrature, the upper limit of 
the integration of Eq. (fT7|) should be assigned. We decided to take z max = 16. 5x + 3.5 so 
that the integrand is small enough at z — z max for the x region 10~ 5 < x < 0.99. If the 
number Nql is small, the evolved nonsinglet distribution is not accurate enough in the 
small and large x regions as shown in Fig. HJ 
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Figure 4. Evolved nonsinglet distribution ra- 
tios x(u v + d v ) Melhn /x(u v + d v ) BF are shown 
for N GL =6, 10, 20, and 50 in the Mellin- 
transformation method. 



Figure 5. Integrand of the Mellin inversion in 
Eq. lO. 



The inaccuracy in the small and large x regions for Ngl=6 and 10 is understood in 
the following way. We show the integrand of Eq. (fTTjl in Fig. El by taking x=10~ 5 , 10~ 3 , 
10 _1 , and 0.9. It is clear that the integrand oscillates at small x so that a certain number 
of Gauss-Legendre points is needed for getting an accurate evolution. On the other hand, 
the integrand does not decrease rapidly at large z for the large-x case (x=0.9), so that 
large z max should be taken for the integration. In addition, the positive contribution at 
z ~ 1 and the negative one at z ~ 2.5 almost cancel each other, which is another source 
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of numerical inaccuracy. 

The evolution results of the singlet distribution are shown in Figs. El and [7| for the 
Laguerre and Mellin methods, respectively. We notice in Fig. El that the accuracy of the 
singlet evolution is much better than the nonsinglet one in the Laguerre method. It comes 
simply from the functional difference between the nonsinglet and singlet distributions. The 
Laguerre method can be used as an accurate evolution method for the singlet evolution. 
The singlet evolution accuracy for the Mellin method is similar to its nonsinglet ones. If 
the point number Nql is large enough, the evolution becomes accurate. 
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Figure 6. Evolved singlet distribution ratios 

xq Laguerre/ xq BF &re shown for jV iaff = 5, 10, 20, 

and 30 in the Laguerre-polynomial method. 



Figure 7. Evolved singlet distribution ratios 

xq MeUinJ xq BF are shown for N GL =&, 10, 20, 

and 50 in the Mellin-transformation method. 



The gluon evolution results are shown in Figs. |H1 and El for the Laguerre and Mellin 
methods. The Laguerre evolution becomes much more accurate than its singlet-quark 
evolution; however, the Mellin evolution becomes slightly inaccurate at large x. 
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Figure 8. Evolved gluon distribution ratios 

xg Laguerre/ xg BF &re ghown for N Lag = b, 10, 20, 

and 30 in the Laguerre-polynomial method. 



Figure 9. Evolved gluon distribution ratios 

xg MelUnj xg BF are shown for A^ GjL =6, 10, 20, 

and 50 in the Mellin-transformation method. 



4.2. Fragmentation functions 

The fragmentation functions (FFs) are essential for understanding hadron productions 
in high-energy reactions. In addition, they are important for describing ultrahigh-energy 
cosmic rays in the top-down scenario [EJ El- I n comparison with the situation of the 
PDFs, their determination is still premature in the sense that experimental data are not 
still enough to determine them accurately. However, there are available e + e~ annihilation 
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data, which could be used for a global analysis of the FFs. The current status of such 
analyses is summarized in Ref. [|2]. We use the KKP parametrization [E3 at Q 2 =2 GeV 2 
as the initial functions, and then they are evolved to Q 2 =100 GeV 2 with the KKP scale 
parameter for testing three evolution methods. 

The brute-force evolution with N t =200 and N x =A000 is taken as the standard for 
showing other evolution results as it was done in the previous subsection. In Fig. 
the Laguerre-method results for the singlet fragmentation function into the proton and 
antiproton, = Y^A^q^ + -^fe )> are shown. The Mellin method results are shown 

in Fig. ^2 The small- £ region is shown in these figures for comparing the results with 
the PDF accuracy in Sec. 14. 1\ although it is outside the range of current accelerator 
experiments. As it was found in the PDF evolution, the Laguerre method is not excellent 
in the small- and large-x regions unless a large number of polynomials is taken. The ratios 
in Figs. El and ^2 show a similar tendency to the ratios of the PDF singlet evolution 
results in Figs. Eland[7[ respectively. 
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The gluon FF evolution results are shown in Figs. and EH for the Laguerre and 
Mellin methods. The gluon FF evolution by the Laguerre method is accurate in the 
most-x region except for the large- a; part. The Mellin method is also accurate except for 
the large- a; region. 

We have compared three evolution methods; however, there are other methods [ I16j. 
Because the numerical solution of the DGLAP equations is very important for describing 
high-energy hadron reactions, an efficient and accurate method should be investigated 
further. 

4.3. Computation time 

We show typical computation time for each evolution method. However, we should 
aware that it depends much on the numbers, N t and N x in the brute-force method, Ni ag 
in the Laguerre method, and Nql in the Mellin method. Therefore, we list CPU time for 
different parameter values for N t , N x , N Lag , and N GL in Tabled] by running the codes for 
the nonsinglet PDF. In the singlet evolution, the time becomes longer but the tendency of 
three evolution methods is the same. The used machine is DELL-Dimension-8800 with a 
Pentiam-4 2.8G CPU. The operating system is Redhat-Linux 8.0 and the fortran compiler 
is g77. 

Table 1 

CPU time for calculating PDFs at 500 x-points by running each evolution code for the 
nonsinglet distribution with the linux-g77 compiler on a Pentiam-4 machine with a 2.8G 
CPU. 



Method 


Condition 


CPU time (seconds) 
for PDFs at 500 appoints 


Brute-force 


JV t =50, N x =1000 


1.501 


N t =200, N x =1000 


5.986 


iV t =200, N x =4000 


95.634 


Laguerre 


N Lag =5 


0.005 


N La9 =10 


0.011 


N La9 =20 


0.025 


N Lag =30 


0.044 


Mellin 


N GL =6 


0.154 


N GL =10 


0.244 


N GL =20 


0.464 


N GL =50 


1.128 



Among the three methods, the brute-force method takes the longest time for the compu- 
tation simply because the large numbers of steps are taken. Therefore, it typically takes a 
few seconds for obtaining a reasonable accuracy for the nonsinglet evolution (A^=50-200, 
A^=1000). 

In the Laguerre method, the evolution is simply given by the summation of the La- 
guerre coefficients which are calculated partially with the recursion relation. There is no 
numerical integration involved in the evolution calculation, so that this method is by far 
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the fastest among the studied methods. Even if N Lag =30 is taken, it takes 0.044 seconds 
for the nonsinglet evolution. It means that it is one hundred times faster than the brute- 
force method. If one is interested in using it for the singlet evolution and if one does not 
mind one or two percent error, it is certainly the best method. 

In the Mellin method, accurate evolution results are obtained with Nql — 20. The 
computation time is significantly shorter and it is several times faster than the brute- 
force method. This is the reason why this method is popular among high-energy physics 
researchers. 

5. Summary 

We have compared the evolution results of the parton distribution functions and the 
fragmentation functions by using three evolution methods, brute-force, Laguerre-polynomial, 
Mellin-transformation methods. The advantages and disadvantages of each method are 
summarized in Tabled 



Table 2 

Summary of advantages and disadvantages of each evolution method. 



Method 


Advantage 


Disadvantage 


Brute-force 


Simple code: The computer code is 
very simple. More complicated evo- 
lution equations with higher-twists 
{e.g. in Ref. [HI]) could be han- 
dled easily. The evolution could be 
accurate in the small- and large- x re- 
gions. 


Long computation time: In order to 
obtain an accurate evolution, large 
numbers of steps {N t and N x ) are 
needed. If one uses the code for 
many evolution calculations, it takes 
a significant amount of time. 


Laguerre 


Very fast: It takes less than a second 
by an ordinary desktop computer. 
As long as one does not mind the 
very small- and large- x regions, it is 
a good method for repeated evolu- 
tion calculations. 


Accuracy at small and large x: De- 
pending on the initial functional 
form, the results do not converge 
at small x unless a large number of 
polynomials are used. It is also diffi- 
cult to obtain accurate evolution at 
large x. 


Mellin 


Fast: By choosing an appropriate 
z max and Nql in each x region, the 
code becomes much faster than the 
brute-force computation. For re- 
peated evolution calculations with 
certain accuracy, this method is ap- 
propriate. 


Accuracy at small and large x: One 
should be careful about the choices 
of z max and Nql- In particular, the 
Mellin inversion should be carefully 
done at very large x. 
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